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Abstract 



The use of reaction-diffusion models rests on the key assumption that the 
underlying diffusive process is Gaussian. However, a growing number of stud- 
ies have pointed out the presence of anomalous diffusion, and there is a need 
to understand the dynamics of reactive systems in the presence of this type 
of non-Gaussian diffusion. Here we present a study of front dynamics in 
reaction-diffusion systems where anomalous diffusion is due to asymmetric 
Levy flights. Our approach consists of replacing the Laplacian diffusion op- 
erator by a fractional diffusion operator, whose fundamental solutions are 
Levy a-stable distributions. Numerical simulation of the fractional Fisher- 
Kolmogorov equation and analytical arguments show that anomalous diffu- 
sion leads to the exponential acceleration of fronts and a universal power law 
decay, x~ a , of the tail, where a, the index of the Levy distribution, is the 
order of the fractional derivative. 
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Reaction-diffusion models have found widespread applicability in many areas, including 
chemistry, biology, physics and engineering; see, for example, Refs. [1] and references therein. 
The simplest reaction-diffusion models are of the form 

dt<j> = x%<t> + F(<i>), (1) 

where x is the diffusion constant, and F is a nonlinear function representing the reaction 
kinetics. Examples of particular interest include the Fisher-Kolmogorov equation for which 
F = 70(1 — 0) and the real Ginzburg-Landau equation for which F = 70(1 — 2 ). The non- 
trivial dynamics of these systems arises from the competition between the reaction kinetics 
and diffusion. 

At a microscopic level, diffusion is the result of the random motion of individual particles, 
and the use of Laplacian operators, <9^0, to model it rests on the key assumption that this 
random motion is an stochastic Gaussian process. However, a growing number of works have 
shown the presence of anomalous diffusion processes for which the mean square variance, 
([x — (x)] 2 } ~ t 7 , grows faster (7 > 1, in the case of superdiff usion) , or slower (7 < 1, in the 
case of subdiffusion) than in a Gaussian diffusion process [2] . Accordingly, an important open 
problem is to understand the dynamics of reaction-diffusion systems when the assumption 
of Gaussian diffusion fails. This problem has a particular relevance to plasma physics, 
which was our original motivation to carry out this study. In particular, reaction-diffusion 
models with Gaussian diffusion have been used to study the dynamics of spatio-temporal 
propagating fronts in the transition to high confinement regimes in magnetically confined 
plasmas [3]. However, these studies must incorporate the fact that, in this type of plasmas, 
evidence of non-Gaussian diffusion has been observed in perturbative transport experiments 
[4], in numerical simulation of three dimensional turbulence [5], and in test-particle transport 
studies [6]. 

The origin of non-Gaussian diffusion can be traced back to the existence of long-range 
correlations in the dynamics, or the presence of anomalously large particle displacements 
described by broad probability distributions. Here we are interested in the second possibility. 
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In particular, we focus on systems that exhibit anomalous diffusion caused by Levy flights, 
for which the probability distribution of particle displacements, p(£), is broad in the sense 
that (£ 2 ) = oo. As it is well-known, for these kind of systems the central limit theorem 
cannot be applied; and as iV — > oo, the probability distribution function of x = J2n 
rather than being Gaussian, is an a-stable Levy distribution [2,7] 

The majority of studies on anomalous diffusion due to Levy flights have focused on sym- 
metric processes for which p(£) = p(—£). However, this is not always the case. For example, 
numerical studies of test particle transport in magnetized plasmas and in geophysical flows 
show that these systems have a built-in asymmetry that gives rise to asymmetric transport 
[6,8]. Also, non-Gaussian asymmetric processes are likely to have an important role in the 
asymmetries observed in pulse propagation experiments in confined plasmas [4]. Motivated 
by this, we focus on anomalous diffusion processes that exhibit Levy flights in one direction, 
say for x > 0, but Gaussian behavior in the other direction, x < 0. The main goal of this 
paper is to understand how this asymmetry manifests in the propagation of right-moving 
(x > 0) and left-moving (x < 0) fronts in reaction-diffusion systems. To study this problem 
we propose the following model: 



where the Laplacian operator d 2 has been replaced by a D" , the Riemann-Liouville, left frac- 
tional derivative of order a, with 1 < a < 2. Fractional calculus is a natural mathematical 
generalization of standard calculus [9], for example, -^D" e tkx = (ik) a e tkx , and in recent 
years has been applied to a large class of problems in the applied sciences; see Refs. [10,11] 
and references therein. 

Previous studies of reaction-diffusion systems with anomalous diffusion include Ref. [12], 
where reaction in the presence of subdiffusion was studied using fractional derivative oper- 
ators in time. In contrast, here we are interested in superdiffusion, which we model with 



dt<l> = XaDl4> + F{4>), 



(2) 




(3) 
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fractional operators in space. The interplay of bistable reaction processes and anomalous 
diffusion caused by Levy flights was addressed in Ref. [13]. More recently, superfast front 
propagation in reactive systems with non-Gaussian diffusion was discussed in Ref. [14] . The 
model studied in Ref. [14] consisted of a time-discrete reaction system coupled to a superdif- 
fusive Levy process described by an integral operator with an algebraic decaying propagator. 
In our model in Eq. (2), we do not assume a time-discrete reaction kinetics; and through 
the use of fractional operators we use the exact Levy propagator. In addition, we consider 
here the role of asymmetric transport. 

The physical motivation for using the left-fractional derivative rests on the fact that for 
F — 0, Eq. (2) reduces to an asymmetric, space-fractional diffusion equation whose solution 
for a delta function initial condition, <f>(x, t = 0) = S(x), in the infinity domain x G (— oo, oo) 
is 

where 

PM = ^ r ^^dk. (5) 

For a = 2, Eq. (5) is the Gaussian propagator, and Eq. (4) is the fundamental solution of 
the standard diffusion equation. However, for 1 < a < 2, Eq. (5) is an extremal, a-stable 
Levy distribution, i.e., a distribution with maximum skewness. These distributions are the 
attractors of stochastic processes that exhibit Levy flights in only one direction, they have 
algebraic asymptotic behavior at plus infinity, p a (i]) ~ l/r) a+1 , but decay exponentially at 
minus infinity [7]. In principle, one could add to Eq. (2) a right-fractional derivative with 
a suitable weighting factor, and for F = obtain Levy a-stable distribution of arbitrary 
skewness [15]. Equation (4) implies that (x n ) = (xt) n ^ a f-rV n Pa(v)dv diverges as r — > oo, 
for n > a. However, in physical applications (e.g., Ref. [8]), a finite-r cut-off leads to the 
finite-size scaling (x n ) ~ t n ^ a , which implies super diffusive behavior with 7 = 2/ a. 

In the remainder of this paper, we present a numerical and analytical study of front 
solutions of Eq. (2) with a reaction dynamics of the Fisher-Kolmogorov type, F = 70(1 — 0), 
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(xt)v« 



(4) 



with two fixed points, one stable, = 1, and the other unstable, = 0. As it is usually done 
for initial value problems with fractional derivatives in time, and boundary value problems in 
finite domains with fractional derivatives in space [10,9], we regularize the singular behavior 
of the fractional operator by subtracting the value of <f>(x) at the lower limit, namely 



where 0(a) = <f>(x = a,t). In an infinite domain, a — > — oo, and Eq. (6) reduces to Eq. (2) 
since the fractional derivative, -oo-D", of a constant is zero. Here we consider the finite 
domain, x G (0, 1), and take a = 0. The numerical solutions were obtained by integrating 
Eq. (6) with boundary conditions 0(0) = 1 and 0'(O) = for right-propagating fronts, and 
0'(O) = and 0(1) = 1 for left-propagating fronts. We used a semi-implicit time advance 
with an up-wind finite difference scheme. The fractional operator was discretized using 
the Grunwald-Letnikov definition of the fractional derivative [9]. Details on the numerical 
method will be published elsewhere. 

We consider two classes of initial conditions 



where 4>q(x) takes the — sign and corresponds to a right-propagating front, and <f> l (x) takes 
the + sign and corresponds to a left-propagating front. Figure 1 shows the time evolution 
of the front profile for 4> r (x) with W = 0.001 and x = 0.003, obtained from the direct 
numerical integration of Eq. (6) with 7 = 1 and x = 5 x 10~ 7 . In this case, the front 
propagates to the right and develops an algebraic decaying tail ~ x~ a . The time evolution 
of the tail, <p(x = l,t), exhibits, for large t, exponential growth ~ e 7 *. The acceleration of 
the front is evident in the space-time diagram in Fig. 2 that shows a contour plot of <f>(x,t). 
Front acceleration and algebraic decay of the tail was also observed in the model studied in 
Ref. [14], but with a different exponent, namely ~ x~( a+1 \ 

Left-propagating fronts exhibit a more standard dynamics. In particular, Fig. 3 shows 
the time evolution of the front profile for an initial condition (f) l (x) with W = 0.001, 7 = 1, 



dt<f> = XaD^ [0- 0(a)] +70 (1-0), 



(6) 




(7) 
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and x — 5 x 1CT 7 (the same parameter values used in Fig. 1) and x = 0.9. In this case, the 
front exhibits an exponential decay and a self-similar propagation with constant speed c. 

The numerical results presented above can be explained analytically using the leading- 
edge approximation, extensively used in the study of fronts with Gaussian diffusion [16]. In 
this approximation, in the limit 1, the reaction kinetics F(cf>) is linearized around the 
unstable phase leading to the linear fractional equation: 

d t <f> = X -ooD^(f> + j(f>. (8) 

Substituting = e yt tp(x,t) into Eq. (8) yields a fractional diffusion equation for -0 whose 
general solution is 

if>(x,t)= / Pa ( V )Mx-(xt) 1/a v}dv, (9) 



where, as discussed before, p a (rj) given in Eq. (5) is the extremal Levy a-stable distribution 
that is the Green's function of the asymmetric fractional diffusion equation. For a front 
initial condition of the form (f) (x < 0) = 1 and (f) (x > 0) = e~ Xx , Eq. (9) gives 

rx( X tr 1/a 
lx(xt)- x / a J-oo 

Before using this solution to obtain the asymptotic behavior of fractional diffusion, right- 
propagating fronts, it is instructive to consider the standard diffusion limit, a — 2. In this 
case, the Green's function is the Gaussian propagator P2(v), and Eq. (10) becomes 



(/>(x,t) = e 7 * / Pa (r]) dr] + e~ Xx+ ^ / e x{xt)1/a " p a ( v ) drj . (10) 



1-P< X 



(11) 



where P(z) is the normal probability distribution function. Using the fact that P(z — > oo) = 
1, we get from Eq. (11) the asymptotic behavior <j>(x,t) ~ e -A(x-ct)_ where the speed of 
the front, c, is related to the steepness of the front, A, according to c = 7/A + %A, with the 
minimum front speed c m = selected for sufficiently steep profiles [1,16]. 

In the fractional case we consider a large, fixed t, and x(xt)~ 1 ^ a — > 00. Introducing a 
cut-off Q such that 1 Q < x(%t) _1 / a in the second integral in Eq. (10), integrating by 
parts, and using the asymptotic expression of the Levy distribution, p a (v) ~ l/i] a+1 we get 



Xte 



it 



a 



(1 + a) r<xtr 1/a e A M v ^H 

+ — + 



tjC 



7^ L ^ ^ + • • • 



(12) 



A \( x t) 1+1/a Jn r? 2 

where the dots denote terms of order e~ Xx . Since the integrand on the right hand side of 
Eq. (12) is bounded by l/r] 2+a , the third term in this equation is at most of order x _a_1 . 
Therefore, in agreement with the numerical result in Fig. 1, to leading order the tail of right 
propagating fronts decays as ~ l/x a . The time-asymptotic dynamics for fixed, large x 
can also be obtained from Eq. (10). In this case, a stationary phase approximation gives, in 
agreement with the numerical results, ~ e 7 * to leading order, with a correction of order 

The previously discussed asymptotic results can be summarized as ~ x~ a e 7 *, which, 
as Fig. 2 shows, reproduces the numerical results. Solving the equation 0(x, £) = 0^ for 
x, for a fixed value = 0l G (0,1) gives the Lagrangian trajectory xl = x(t; 0l) of the 
coordinate of a point in the front with concentration = <p L . In the asymptotic limit we 
have xl ~ e 7t//Q . That is, the Lagrangian velocity vl = dxi/dt of right-moving fronts grows 
exponentially with time, vl ~ e 7t//a . 

To conclude, consider the dynamics of left-moving fronts, like the one shown in Fig. 3. 
In this case, the leading-edge description is straightforward since these fronts exhibit expo- 
nential tails and propagate at constant speed. Substituting ~ exp[A(rc + ct)} into Eq. (2), 
we obtain the dispersion relation c = c(A), with minimal front speed c m = c(A m ) where 

i , x „ _ ./« I_2_ ^ ( °- 1)/<> 



i + c m = a X ^ [ -^- l) . (13) 

Numerical results support the idea that for steep (A > A m ) initial conditions the front 
selects the minimum velocity c m , whereas for wide (A < A m ) initial conditions the front 
speed depends on the initial condition according the first equation in (13). As expected, 
in the limit a = 2 the well-known results of front propagation in the presence of Gaussian 
diffusion [1,16] are recovered. 

Summarizing, in this paper we have proposed the use of fractional-diffusion operators 
to study front dynamics in reaction-diffusion systems with non-Gaussian diffusion caused 
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by asymmetric Levy flights. Numerical and analytical results show that right-moving fronts 
accelerate exponentially, and develop an algebraic decaying tail, <fi ~ x~ a e lt . Left-moving 
fronts have exponential decaying tails and move at a constant speed given by Eq. (13). The 
results are general in the sense that they are independent of the details of the reaction ki- 
netics, provided it is of the "pull" type with a stable and an unstable phase. Two areas 
of potential application of the ideas presented here are Plasma Physics and Biology. In 
particular, transport studies in three-dimensional pressure-driven plasma turbulence [5] and 
in drift waves [6] have shown evidence of anomalous diffusion and non-Gaussian probability 
distributions of particle displacements which can be modeled using fractional diffusion equa- 
tions. Also, fractional diffusion equations seem to capture important aspects of perturbative 
transport experiments in fusion plasmas including nonlocal diffusion effects in the propa- 
gation of cold pulses (e.g., [4]). On a parallel development, reaction-diffusion models with 
Gaussian diffusion have been used to study the turbulence- shear flow interaction in the L-H 
transition in fusion plasmas (e.g., [3]), and there is a pressing need to understand the role of 
non-Gaussian, non-local diffusion in the dynamics of L-H transition fronts. The results pre- 
sented here represent a first step in the study of this open problem. On the other hand, the 
fractional Fisher-Kolmogorov equation discussed here shares some similarities with reaction 
equations with integro-differential operators used in Biology to model long-range diffusion 
and spatial patterns in neural firing [1]. In this regard, the algebraic decaying kernel in the 
fractional diffusion operator might be useful for describing strongly non-local processes. 

We thank Angelo Vulpiani and Michael Menzinger for useful conversations. This work 
was sponsored by the Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for 
the U.S. Department of Energy under contract DE-AC05-00OR22725. 



9 



REFERENCES 

[1] J. D. Murray, Mathematical Biology (Springer Verlag, New York, 1989); M. C. Cross 
and P. C. Hohenberg, Rev. Mod. Phys., 65, 851 (1993). 

[2] J. P. Bouchaud and A. Georges, Phys. Rep. 195, 127 (1990). 

[3] D. del-Castillo-Negrete, and B. A. Carreras, Phys. Plasmas, 9, 118 (2002); D. del- 
Castillo-Negrete, B. A. Carreras, and V. Lynch, Physica D, 168-169, 45-60 (2002). 

[4] K. W. Gentle, et.al., Phys. Plasmas, 2, 2292 (1995). 

[5] B. A. Carreras, V. E. Lynch, and G. M. Zaslavsky, Phys. Plasmas 8, 5096 (2001). 
[6] D. del-Castillo-Negrete, Phys. Plasmas 7, 1702 (2000). 

[7] G. Samorodnitsky, and M. S. Taqqu, Stable non-Gaussian random processes (Chapman 
& Hall, New York, 1994). 

[8] T. H. Solomon, E. R. Weeks, and H. L. Swinney, Phys. Rev. Lett. 71, 3975 (1993); D. 
del-Castillo-Negrete, Phys. Fluids 10, 576 (1998). 

[9] I. Podlubny, Fractional Differential Equations (Academic Press, San Diego, 1999). 

[10] R. Metzler, and J. Klafter, Phys. Rep., 339, 1, (2000). 

[11] R. Hilfer, Application of Fractional Calculus in Physics, (World Scientific, Singapore, 
2000); I. M. Sokolov, J. Klafter, and A. Blumen, Phys. Today 55, 48 (2002). 

[12] B. I. Henry, and S. L. Wearne, Phys. A, 276, 448, (2000). 

[13] D. H. Zanette, Phys. Rev. E, 55, 1181, (1997). 

[14] R. Mancinelli, D. Vergni, and A. Vulpiani, Europhys. Lett., 60, 532-538, (2002). 

[15] A. I. Saichev, and G. M. Zaslavsky, CHAOS, 7, 753, (1997); A. S. Chaves, Phys. Lett. 
A, 239, 13-16 (1998); P. Paradisi, R. Cesari, F. Mainardi, and F. Tampieri, Physica A, 
293, 130-142 (2001). 

10 



[16] W. van Saarloos, Phys. Rev. A 37, 211 (1988). 



11 



FIGURES 
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FIG. 1. Right propagating front profiles at successive times, obtained from a numerical integra- 
tion of the fractional Fisher-Kolmogorov Eq. (6) with a = 1.5 and initial condition <f>(x, 0) = <^>o(x) 
in Eq. (7). The dashed line has slope equal to a. In agreement with the analytical result, 
right-moving fronts develop an algebraically decaying tail, <p ~ x~ a - 
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FIG. 2. Contour plot of the numerical solution (f>(x,t) in Fig. 1. The curvature of the 
iso-contours illustrates the exponential acceleration of fronts in the fractional Fisher-Kolmogorov 
equation. The dashed line corresponds to the analytical scaling result <p ~ x~ a e lt . 
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FIG. 3. Left-propagating front profiles at successive times obtained from a numerical integra- 
tion of the fractional Fisher-Kolmogorov Eq. (6) with a = 1.5 and initial condition <j>(x, 0) = (j> l {x) 
in Eq. (7). The dashed line has slope equal to A = 2/W = 2000. In agreement with the analytical 
result in Eq. (13), the front exhibits an exponential decay and has a speed of c = 5.22 x 10~ 4 . 
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